/**********************************************************
 *
 *  Copyright (c) 2003  SeikoEpson Inc
 *  All Rights Reserved
 *
 *  File name : f_atan2.c
 *  Function  :
 *        C unStandard Library Function "f_atan2"
 *        The "f_atan2" function calculates the arc tangent of "y/x".
 *
 *        This file is copied from math.lib of 
 *                  CC33 tool(CC33v40.exe + CCPA45.exe, 
 *                  math of ansilib33v25 )
 *  original file's Revision  :
 *      1999/12/08    Created by H.Ishida
 *      2000/01/06    Fixed a bug by H.Ishida
 *      2000/02/14    change value M.Igarashi
 *  Revision  :
 *      2003/04/08    IrumaSoft M.Takeishi   1.st design
 *
 **********************************************************/
/* -- Format of floating point (single precision) ---
 * 
 * 31|30         23|22                0|
 *  -----------------------------------
 * | |   exponent  |   fraction        |
 *  -----------------------------------
 * |    8 bits          23 bits
 */


#include <f_math.h>
#include <f_smcvals.h>

#define	K0		0.445413400592906803197511e2f
#define	K1		0.77477687719204208616481e2f
#define	K2		0.40969264832102256374186e2f
#define	K3		0.666057901700926265753e1f
#define	K4		0.1589740288482307048e0f
#define	S0		0.445413400592906804445995e2f
#define	S1		0.92324801072300974840693e2f
#define	S2		0.62835930511032376833267e2f
#define	S3		0.1550397755142198752523e2f

float f__atan( float );

/******************************************************************************
 *	Function: float atan2(float sfX, float sfY)
 ******************************************************************************/
float f_atan2(float sfX, float sfY)
{
	long lNanData;
	long lX,lX2;
	float	sfTemp;
	float sfX2,sfRet;
	
	F_GETW(lX,sfX);				// get x
	F_GETW(lX2,sfY);			// get y



	/* Is NaN about sfX? */
	lNanData = 0x7fc00000;
	if ( ( lX & lNanData ) == lNanData || 
	     ( lX2 & lNanData ) == lNanData ) {
	     errno = EDOM;
	     sfRet = f_NAN._F;
		return sfRet;
	}

	sfTemp=sfX+sfY;
	if (sfTemp==sfX){
		
		if ( ( lX == 0x00000000 ) && ( lX2 == 0x00000000 ) ) {
			/* if x = y = 0.0f, errno should be EDOM. */
			errno = EDOM;
		}

		sfRet=F_HLFP;
		if( ( lX & F_NEG_ZERO_L ) == 0 ) {
			return sfRet;
		} else {
			return -sfRet;
		}
	}else if( (lX2&0x80000000) != 0 ) {

		F_SET_SIGN( sfX2, sfX );	// clear sign

		sfRet = F_WHP - f__atan(sfX2/sfY);
		//if( (lX&0x80000000) == 0 || ( lX ==0x0 && ulHx==0x0 )){
		if( (lX&0x80000000) == 0 ){     // why ( lX ==0x0 && ulHx==0x0 ) exists???
			return sfRet;
		} else {
			return -sfRet;
		}
	}
	//else if((lX&0x80000000) == 0 || ((lX&0x80000000) == 0) && (ulHx !=0x0 )){
	else if( (lX&0x80000000) == 0 ){   // why   ((lX&0x80000000) == 0) && (ulHx !=0x0 )   exists...?
		sfRet = f__atan(sfX/sfY);
		return sfRet;
	} else {
		sfRet = -f__atan(-sfX/sfY);
		return sfRet;
	}
}

/******************************************************************************
 *	Function: float f_atan(float sfX)
 ******************************************************************************/


float f__atan(float sfX)
{
	long lX, lNanData;
	float sfXsq,sfRet,sfTemp;
	

	F_GETW(lX,sfX);		// get low


	/* Is NaN ? */
	lNanData = 0x7fc00000;
	if ( (lNanData & lX ) == lNanData ) {
		errno = EDOM;
		sfRet = f_NAN._F;
		return sfRet;
	}
	
	if( (0x80000000 & lX ) != 0 || lX < 0x3ed413cd ){	// check qarea  0x3ed413cd = 0.412....
		sfXsq = sfX*sfX;
		sfRet=K4*sfXsq+K3;
		sfRet=sfRet*sfXsq+K2;
		sfRet=sfRet*sfXsq+K1;
		sfRet=sfRet*sfXsq+K0;
		sfTemp=sfXsq+S3;
		sfTemp=sfTemp*sfXsq+S2;
		sfTemp=sfTemp*sfXsq+S1;
		sfTemp=sfTemp*sfXsq+S0;
		sfRet = sfRet/sfTemp;
		sfRet = sfRet * sfX;
		return sfRet;
	}else if( (0x80000000 & lX ) == 0 && ( lX > 0x409a827a )){	// check qarea   0x409a827a = 4.75....
	
		sfXsq = 1.0f / sfX;
		sfXsq = sfXsq*sfXsq;
		sfRet=K4*sfXsq+K3;
		sfRet=sfRet*sfXsq+K2;
		sfRet=sfRet*sfXsq+K1;
		sfRet=sfRet*sfXsq+K0;
		sfTemp=sfXsq+S3;
		sfTemp=sfTemp*sfXsq+S2;
		sfTemp=sfTemp*sfXsq+S1;
		sfTemp=sfTemp*sfXsq+S0;
		sfRet = F_HLFP-sfRet/sfTemp*(1.0f/sfX);
		return sfRet;


	}else{

		sfXsq=(sfX-1.0f)/(sfX+1.0f);
		sfXsq=sfXsq*sfXsq;
		sfRet=K4*sfXsq+K3;
		sfRet=sfRet*sfXsq+K2;
		sfRet=sfRet*sfXsq+K1;
		sfRet=sfRet*sfXsq+K0;
		sfTemp=sfXsq+S3;
		sfTemp=sfTemp*sfXsq+S2;
		sfTemp=sfTemp*sfXsq+S1;
		sfTemp=sfTemp*sfXsq+S0;
		sfRet = F_QUAP+sfRet/sfTemp*((sfX - 1.0f) / (sfX + 1.0f));
		return sfRet;

	}


}
